#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>

/* Default number of iterations */
#define ITERATIONS 25

/* Default minimum precision, in bits */
#define PRECISION 35000000

int
main(int argc, char *argv[])
{
    int i, iterations, precision;
    char *p_str;
    mp_exp_t exp;
    mpf_t x, x0, y, y0, p, p0;
    mpf_t tmp, sqrt_x0, inv_sqrt_x0;

    iterations = ITERATIONS;
    /* The number of iterations; if given but is not invalid, abort */
    if (argc > 1 && sscanf(argv[1], "%d", &iterations) != 1)
        exit(1);

    precision = PRECISION;
    /* The precision; if given but is not valid, abort */
    if (argc > 2 && sscanf(argv[2], "%d", &precision) != 1)
        exit(1);

    mpf_set_default_prec(precision);
    mpf_inits(x, x0, y, y0, p, p0, NULL);
    mpf_inits(tmp, sqrt_x0, inv_sqrt_x0, NULL);

    mpf_set_d(x0, 2.0);
    mpf_sqrt(x0, x0); /* x0 = sqrt(2) */
    mpf_sqrt(y, x0); /* y = sqrt(sqrt(2)) */
    mpf_add_ui(p0, x0, 2); /* p0 = sqr(2) + 2 */

    /* Calculating x from x0 */
    mpf_sqrt(sqrt_x0, x0); /* sqrt_x0 = sqrt(x0) */
    mpf_ui_div(inv_sqrt_x0, 1, sqrt_x0); /* inv_sqrt_x0 = 1 / sqrt_x0 */
    mpf_add(x0, sqrt_x0, inv_sqrt_x0); /* x0 = sqrt_x0 + inv_sqrt_x0 */
    mpf_div_ui(x, x0, 2); /* x = x0 / 2 */

    /* Store x into x0 for the next iteration */
    mpf_set(x0, x);

    /* Calculating p from p0, x and y */
    mpf_set(y0, y);
    mpf_add_ui(x, x, 1);
    mpf_add_ui(y, y, 1);
    mpf_div(tmp, x, y);
    mpf_mul(p, p0, tmp);

    /* Store p into p0 for the next iteration */
    mpf_set(p0, p);

    for (i = 1; i < iterations; ++i) {
        /* Calculating x from x0 */
        mpf_sqrt(sqrt_x0, x0); /* sqrt_x0 = sqrt(x0) */
        mpf_ui_div(inv_sqrt_x0, 1, sqrt_x0); /* inv_sqrt_x0 = 1 / sqrt_x0 */
        mpf_add(x0, sqrt_x0, inv_sqrt_x0); /* x0 = sqrt_x0 + inv_sqrt_x0 */
        mpf_div_ui(x, x0, 2); /* x = x0 / 2 */

        /* Store x into x0 for the next iteration */
        mpf_set(x0, x);

        /* Calculating y from y0 */
        mpf_mul(tmp, y0, sqrt_x0); /* tmp = y0 * sqrt_x0 */
        mpf_add(tmp, tmp, inv_sqrt_x0); /* tmp = tmp + inv_sqrt_x0 */
        mpf_add_ui(y0, y0, 1); /* y0 = y0 + 1 */
        mpf_div(y, tmp, y0); /* y = tmp / y0 */

        /* Store y into y0 for the next iteration */
        mpf_set(y0, y);

        /* Calculating p from p0, x and y */
        mpf_add_ui(x, x, 1);
        mpf_add_ui(y, y, 1);
        mpf_div(tmp, x, y);
        mpf_mul(p, p0, tmp);

        /* Store p into p0 for the next iteration */
        mpf_set(p0, p);
    }

    p_str = mpf_get_str(NULL, &exp, 10, 0, p);
    printf("%s\n", p_str + exp);
    free(p_str);

    mpf_clears(tmp, sqrt_x0, inv_sqrt_x0, NULL);
    mpf_clears(x, x0, y, y0, p, p0, NULL);

    return 0;
}
